--- title: Peak slice maps (working on it) keywords: fastai sidebar: home_sidebar summary: "Slicing the cube " description: "Slicing the cube " nb_path: "notebooks/60_peakmaps.ipynb" ---
import skimage.exposure as ske
from maxrf4u import DataStack, HotmaxAtlas, get_hotslices, get_peakmaps, multi_plot, plot_peakslices
ds = DataStack('RP-T-1898-A-3689.datastack')
x_keVs = ds.read('maxrf_energies')
y_max = ds.read('maxrf_maxspectrum')
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
imvis_reg_highres = ds.read('imvis_reg_highres')
imvis_reg = ds.read('imvis_reg')
imvis_extent = ds.read('imvis_extent')
hma = HotmaxAtlas('RP-T-1898-A-3689.datastack')
hma.hotmax_pixels[:,2]
def hotslice_indices(hma=hma):
'''Temporary hack to locate hotpixel slice index for (indicated with a star in the multiplot).
In the future do better ordering by putting hotslice first.'''
keV_idxs = hma.hotmax_pixels[:,2]
hotslice_idxs = []
for n, peak_idxs in enumerate(hma.peak_idxs_list):
hotslice_idxs.append(peak_idxs.index(keV_idxs[n]))
return hotslice_idxs
hotslice_indices()[9]
def hotslice_index(n, hma=hma):
'''Temporary hack to locate hotpixel slice (indicated with a star).
In the future do better ordering by putting hotslice first.'''
keV_i = hma.hotmax_pixels[:,2][n]
hotslice_i = hma.peak_idxs_list[n].index(keV_i)
return hotslice_i
hotslice_index(3)
n = 5
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
ax, labels = hma.plot_spectrum(n)
plot_peakslices(x_keVs, slices, ax=ax);
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
titles[n]= '*' + titles[n]
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
fig, axs = multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel)
Let's check if peak #5[6] is an escape peak for Ca_Ka by calculating the energy shift....
peak_idxs = hma.peak_idxs_list[5]
x_keVs[peak_idxs[0]] - x_keVs[peak_idxs[6]]
Indeed!!
title = axs.flatten()[6].get_title() + '(escape)'
axs.flatten()[6].set_title(title)
There is much more to be learned from this overview, but before going into this let's look at other elements.
n = 9 # Mn
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
ax, labels = hma.plot_spectrum(n)
plot_peakslices(x_keVs, slices, ax=ax);
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
Mm, these histogram equalized peak maps aren't very informative. Most of the slices are just noise. I can not be sure how manganese and iron are correlated. Need to plot that in a different way.
FeKa_map = peak_maps[2]
MnKa_map = peak_maps[4]
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])
ax.imshow(MnKa_map)
ax.set_title('Mn_Ka')
y, x, z = hot_pixel
ax.scatter([x], [y], c='r')
ax1.imshow(FeKa_map, vmax=2)
ax1.set_title('Fe_Ka')
ax1.scatter([x], [y], c='r');
My first question is, are there any other manganese speckles? And second, do the correlate with iron speckles, as I would expect?
This question needs to wait. Need speckle segmentation functions.
n = 10 # Fe_Ka
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
ax, labels = hma.plot_spectrum(n)
plot_peakslices(x_keVs, slices, ax=ax);
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
n = 3
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
import moseley as mos
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax)
plot_peakslices(x_keVs, slices, ax=ax);
ax1.set_ylim([0, 1.7])
mos.XFluo('Cl', tube_keV=30).plot(ax=ax1, color='b', peak_labels='full')
mos.XFluo('Rh', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
mos.XFluo('K', tube_keV=30).plot(ax=ax1, color='g', peak_labels='simple')
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
titles[n]= '*' + titles[n]
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
I would like to inspect the hotmax pixel... We need to check if Rh_L3M5 at 2.701 keV is not the major contributor to this
ClKa_map = peak_maps[3]
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])
y, x, z = hot_pixel
ax.imshow(ClKa_map)
ax1.imshow(imvis_reg)
ax.scatter(x, y, c='r')
ax1.scatter(x, y, c='r');
Seems to be a single grain of salt?
n = 15
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
import moseley as mos
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax)
plot_peakslices(x_keVs, slices, ax=ax);
ax1.set_ylim([0, 1.7])
mos.XFluo('Pb', tube_keV=30).plot(ax=ax1, color='k', peak_labels='full')
mos.XFluo('S', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
titles = [f'[{i}]' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
This is most interesting. Let's take a look at the keV map for the overlapping band of sulfur and lead. For sulfur this is the S_KL3 band at 2.311 keV. For lead this is the Pb_M5O2 band and 2.372 keV.
mid_keV = (2.311 + 2.371) / 2
SKa_keV_map = keV_maps[3]
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5], constrained_layout=True)
ax.imshow(peak_maps[3])
ax.set_title('Overlapping band peak_maps[3]')
pos = ax1.imshow(SKa_keV_map, vmin=2.2, vmax=2.45)
fig.colorbar(pos, ax=ax1, shrink=0.8)
ax1.set_title('Energy of maximum \nS=2.311 keV vs Pb=2.372 keV');
Ok, this is the new Gaussian code base. For now we do not return the fancy peak shape parameters to deconvolve overlapping bands. Soon this will be necessary. For now let's simply take a look at the peak maps...